
# Honest

rm(list=ls())

# load packages 
library(foreign)
library(lmtest)
library(sandwich)
library(stargazer)
library(msm)
library("readxl")
library("ggplot2")
library(gridExtra)
library(readstata13)
library(haven)
library(dplyr)
library(forcats)
library(ggeasy)
library(stargazer)
library(tidyverse)
library(broom)
library(ggplot2)
library(RDHonest)
library(dplyr)
library(lmtest)
library(sandwich)
library(RDHonest)
library(rdlocrand)
library(kableExtra)

cep_all <- read_dta("Data/cep_all.dta")


# Table F1
cep_all <- subset(cep_all, 
                  treatment_dosage_coup >= -5 & 
                    treatment_dosage_coup <= 5 & 
                    !is.na(left))

length(cep_all$treatment_dosage_coup[!is.na(cep_all$treatment_dosage_coup)])
length(unique(cep_all$treatment_dosage_coup))
table(cep_all$treatment_dosage_coup)
3340/10

#  Figure F4
barplot <- ggplot(cep_all, aes(x = treatment_dosage_coup)) +
  geom_bar() +
  xlab("Running Variable") +
  ylab("Frequency") +
  theme_minimal()
barplot

ggsave("Output/figuref4.png", plot = barplot, width = 8, height = 6, dpi = 300)

# binomial test

set.seed(7)

# Table F2
binom.test(338, 338+373, 1/2)

library(rdlocrand)

# randomization inference estimation of the minumum window 

# Table F3
bw1 <- rdrandinf(cep_all$left, cep_all$treatment_dosage_coup, wl = -1, wr = 1, seed = 24)


# RD Honest

sample1 <- data.frame()
sample2 <- data.frame()
sample3 <- data.frame()
sample4 <- data.frame()
sample5 <- data.frame()

# Define the range of dosage values
dosage_ranges <- c(1, 2, 3, 4, 5)

# Loop through the dosage ranges and create samples
for (i in dosage_ranges) {
  subset_condition <- cep_all$treatment_dosage_coup >= -i & cep_all$treatment_dosage_coup <= i
  sample_name <- paste("sample", i, sep = "")
  
  # Store the samples in the respective data frame
  assign(sample_name, subset(cep_all, subset_condition))
}

# Table F4: RD Honest (Coefficients and SE were extracted manually)

fit5_rot_nocov <- RDHonest(
  left ~ treatment_dosage_coup,
  data = sample5, cutoff = 0, h = 5,
  kern = "triangular", se.method = "EHW",
  clusterid = cluster)

M_rot <- fit5_rot_nocov$coefficients$M

rd_model5 <- RDHonest(left ~ treatment_dosage_coup | as.factor(encuesta), data = sample5, clusterid = cluster, se.method="EHW", cutoff = 0, h = 5, M = M_rot)
rd_model4 <- RDHonest(left ~ treatment_dosage_coup | as.factor(encuesta), data = sample4, clusterid = cluster, se.method="EHW", cutoff = 0, h = 4, M = M_rot)
rd_model3 <- RDHonest(left ~ treatment_dosage_coup | as.factor(encuesta), data = sample3, clusterid = cluster, se.method="EHW", cutoff = 0, h = 3, M = M_rot)

summ <- function(fit) c(
  est  = unname(fit$coefficients$estimate),
  se   = unname(fit$coefficients$std.error),
  ci_l = fit$coefficients$conf.low, ci_u = fit$coefficients$conf.high, M = fit$coefficients$M,
  N =fit$coefficients$eff.obs)

out <- rbind(`50-60` = summ(rd_model5),
      `51-59` = summ(rd_model4),
      `52-58` = summ(rd_model3))

tab <- out %>%
  as.data.frame() %>%
  rownames_to_column("Window (±h)") %>%
  mutate(
    Estimate   = round(as.numeric(est), 3),
    `Std. Error` = round(as.numeric(se), 3),
    `95% CI`   = sprintf("(%.3f, %.3f)", as.numeric(ci_l), as.numeric(ci_u)),
    M          = round(as.numeric(M), 3),
    `Effective N` = round(as.numeric(N), 3)) %>%
  select(`Window (±h)`, Estimate, `Std. Error`, `95% CI`, M, `Effective N`)

kable(
  tab,
  format   = "latex",
  booktabs = TRUE,
  caption  = "RDHonest estimates (triangular kernel, EHW SEs, survey FEs).") %>%
  kable_styling(latex_options = "hold_position")
